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OO 

Network science is an interdisciplinary endeavor, with methods and appli- 
es 

O cations drawn from across the natural, social, and information sciences. A 

prominent problem in network science is the algorithmic detection of tightly - 
connected groups of nodes known as communities. We developed a general- 
ized framework of network quality functions that allowed us to study the com- 
munity structure of arbitrary multislice networks, which are combinations of 
individual networks coupled through links that connect each node in one net- 
work slice to itself in other slices. This framework allows one to study com- 
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munity structure in a very general setting encompassing networks that evolve 
over time, have multiple types of links (multiplexity), and have multiple scales. 

The study of graphs, or networks, has a long tradition in fields such as sociology and mathemat- 
ics, and it is now ubiquitous in academic and everyday settings. An important tool in network 
analysis is the detection of mesoscopic structures known as communities (or cohesive groups), 
which are defined intuitively as groups of nodes that are more tightly connected to each other 
than they are to the rest of the network (1-3). One way to quantify communities is by a quality 
function that counts intra-community edges compared to what one would expect at random. 
Given the network adjacency matrix A, where the component details a direct connection 
between nodes % and j, one can construct a quality function Q (4, 5) for the partitioning of 
nodes into communities as Q = Y.ij [Aj — Pij] $(9i, 9j), where S(g i: gj) = 1 if the community 
assignments gi and gj of nodes % and j are the same and otherwise, and is the expected 
weight of the edge between % and j under a specified null model. 

The choice of null model is a crucial consideration in studying network community struc- 
ture (2), ideally respecting the type of network studied. After selecting a null model appropriate 
to the network and application at hand, one can use a variety of computational heuristics to 
assign nodes to communities to optimize the quality Q (2, 3). However, such null models have 
not been available for time-dependent networks — one has instead had to use ad hoc methods to 
piece together the structures obtained at different times (6-9) or abandon quality functions for 
an alternative such as the Minimum Description Length principle (10). While tensor decom- 
positions (11) have been used to cluster network data with different types of connections, no 
quality-function method has been developed for such multiplex networks. 

We developed a methodology to remove these limits, generalizing the determination of 
community structure via quality functions to multislice networks that are defined by coupling 
multiple adjacency matrices (see Fig. 1). The connections encoded by the network slices are 
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flexible — they can represent variations across time, across different types of connections, or 
even community detection of the same network at different scales. However, the usual proce- 
dure for establishing a quality function as a direct count of the intra-community edge weight 
minus that expected at random fails to provide any contribution from these inter-slice couplings. 
Because they are specified by common identifications of nodes across slices, inter-slice cou- 
plings are either present or absent by definition, so when they do fall inside communities, their 
contribution in the count of intra-community edges exactly cancels that expected at random. 
In contrast, by formulating a null model in terms of stability of communities under Laplacian 
dynamics, we have derived a principled generalization of community detection to multislice 
networks, with a single parameter controlling the inter-slice correspondence of communities. 

Important to our method is the equivalence between the modularity quality function (12) 
[with a resolution parameter (5)] and stability of communities under Laplacian dynamics (13), 
which we have generalized to recover the null models for bipartite, directed, and signed net- 
works (14). First, we obtained the resolution-parameter generalization of Barber's null model 
for bipartite networks (15) by requiring the independent joint probability contribution to sta- 
bility in (13) to be conditional on the type of connection necessary to step between two nodes. 
Second, we recovered the standard null model for directed networks (16,17) (again with a res- 
olution parameter) by generalizing the Laplacian dynamics to include motion along different 
kinds of connections — in this case, both with and against the direction of a link. By this gener- 
alization, we similarly recovered a null model for signed networks (18). Third, we interpreted 
the stability under Laplacian dynamics flexibly to permit different spreading weights on the dif- 
ferent types of links, giving multiple resolution parameters to recover a general null model for 
signed networks (19). 

We applied these generalizations to derive null models for multislice networks that extend 
the existing quality-function methodology, including an additional parameter uo to control the 
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coupling between slices. Representing each network slice s by adjacencies A; L j S between nodes 
i and j, with inter-slice couplings Cj rs that connect node j in slice r to itself in slice s (see 
Fig. 1), we have restricted our attention to unipartite, undirected network slices (A ijs = A jis ) 
and couplings (Cj rs = Cj sr ), but we can incorporate additional structure in the slices and cou- 
plings in the same manner as demonstrated for single-slice null models. Notating the strengths 
of each node individually in each slice by kj S = J2i M 3 s and across slices by c js = J2 r Cj sr , 
we define the multislice strength by Kj S = kj S + Cj S . The continuous-time Laplacian dynamics 
given by p is = J2j r (Aij S 5 sr + 5ijCj sr )pj r /Kj r — p is respects the intra-slice nature of A ijs and 
the inter-slice couplings of Cj sr . Using the steady state probability distribution p* r = Kj r / (2/x), 
where 2/x = Y,j r K jr, we obtained the multislice null model in terms of the probability p is \j r 
of sampling node % in slice s conditional on whether the multislice structure allows one to step 
from (j, r) to (i, s), accounting for intra- and inter-slice steps separately as 
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The second term in brackets, which describes the conditional probability of motion between 
two slices, leverages the definition of the Cj sr coupling. That is, the conditional probability of 
stepping from (j, r) to (i, s) along an inter-slice coupling is non-zero if and only if % = j, and it 
is proportional to the probability Cj sr /Kj r of selecting the precise inter-slice link that connects 
to slice s. Subtracting this conditional joint probability from the linear (in time) approximation 
of the exponential describing the Laplacian dynamics, we obtained a multislice generalization 
of modularity (see Supporting Online Material for details): 

Qmultislicc ~Z ^ ' \ I Aij s Ts - ^ I &sr ~\~ $ijCj sr \ (5((^j s , Qjr) j 

l P ijsr { \ Zm s / ) 

where we have utilized reweighting of the conditional probabilities, which allows one to have 
a different resolution 7 S in each slice. We have absorbed the resolution parameter for the inter- 
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slice couplings into the magnitude of the elements of Cj sr , which we suppose for simplicity 
take binary values {0, 00} indicating absence (0) or presence (ui) of inter-slice links. 

Community detection in multislice networks can then proceed using many of the same com- 
putational heuristics that are currently available for single-slice networks [though, as with the 
standard definition of modularity, one must be cautious about the resolution of communities (20) 
and the likelihood of complex quality landscapes that necessitate caution in interpreting results 
on real networks (21)]. We studied examples that have multiple resolutions [Zachary Karate 
Club (22)], vary over time [voting similarities in the U.S. Senate (23)], or are multiplex [the 
"Tastes, Ties, and Time" cohort of university students (24)]. We provide additional details for 
each example in the Supplementary Online Material. 

We performed simultaneous community detection across multiple resolutions (scales) in the 
well-known Zachary Karate Club network, which encodes the friendships between 34 members 
of a 1970s university karate club (22). Keeping the same unweighted adjacency matrix across 
slices (A ijs = Aij for all s), the resolution associated to each slice is dictated by a specified 
sequence of ^ s parameters, which we chose to be the 16 values 7 S = {0.25, 0.5, 0.75, . . . , 4}. 
In Fig. 2, we depict the community assignments obtained for coupling strengths u = {0, 0.1, 1} 
between each neighboring pair of the 16 ordered slices. These results simultaneously probe all 
scales, including the partition of the Karate Club into four communities at the default resolution 
of modularity (3, 25). Additionally, we identified nodes that have an especially strong tendency 
to break off from larger communities (e.g., nodes 24-29 in Fig. 2). 

We also considered roll call voting in the United States Senate across time, from the 1st— 
110th Congresses, covering the years 1789-2008 and including 1884 distinct Senator IDs (26). 
We defined weighted connections between each pair of Senators by a similarity between their 
voting, specified independently for each two-year Congress (23). We studied the multislice 
collection of these 110 networks, with each individual Senator coupled to him/herself when 
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appearing in consecutive Congresses. Multislice community detection uncovered interesting 
details about the continuity of individual and group voting trends over time that are simply not 
captured by the union of the 110 independent partitions of the separate Congresses. Figure 3 
depicts a partition into 9 communities that we obtained using coupling to = 0.5. The Con- 
gresses in which three communities appeared simultaneously are each historically significant: 
The 4th and 5th Congresses were the first with political parties; the 10th and 11th Congresses 
occurred during the political drama of former Vice President Aaron Burr's indictment for trea- 
son; the 14th and 15th Congresses witnessed the beginning of changing group structures in 
the Democratic-Republican party amidst the dying Federalist party (23); the 31st Congress 
included the Compromise of 1850; the 37th Congress occurred during the beginning of the 
American Civil War; the 73rd and 74th Congresses followed the landslide 1932 election amidst 
the Great Depression; and the 85th-88th Congresses brought the major American civil rights 
acts, including the Congressional fights over the Civil Rights Acts of 1957, 1960, and 1964. 

Finally, we also applied multislice community detection to a multiplex network of 1640 col- 
lege students at a northeastern American university (24), including symmetrized connections 
from the first wave of this data representing (1) Facebook friendships, (2) picture friendships, 
(3) roommates, and (4) student "housing group" preferences. Because the different connec- 
tion types are categorical, the natural inter-slice couplings connect an individual in a slice to 
him/herself in each of the other 3 network slices. This coupling between categorical slices thus 
differs from that above that connected only neighboring (ordered) slices. Table [T] indicates the 
numbers of communities and the percentages of individuals assigned to 1,2, 3, or 4 communi- 
ties across the four types of connections for different tu, as a first investigation of the relative 
redundancy across the connection types. 

In summary, our multislice framework makes it possible to study community structure in a 
much broader class of networks than was previously possible. Instead of detecting communities 
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in one static network at a time, our formulation generalizing the Laplacian dynamics approach 
of Ref. (13) permits the simultaneous quality-function study of community structure across 
multiple times, multiple resolution parameter values, and multiple types of links. We used this 
method to demonstrate insights in real- world networks that would have been difficult or impos- 
sible to obtain without the simultaneous consideration of multiple network slices. Although our 
examples included only one kind of variation at a time, our framework applies equally well to 
networks that have multiple such features (e.g., time-dependent multiplex networks), and we 
expect multislice community detection to become a powerful tool for studying such systems. 
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Table 1: Communities in the first wave of the multiplex "Tastes, Ties, and Time" network (24), 
using the default spatial resolution (7 = 1) in each of the four slices of data (Facebook friend- 
ships, picture friendships, roommates, and housing groups) under various couplings u> across 
slices, which changed the number of communities and percentages of individuals assigned on a 
per slice basis to 1,2, 3, or 4 communities. 
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Fig. 1: Schematic of a multislice network. Four slices s = {1, 2, 3, 4} represented by adjacen- 
cies Aij S encode intra-slice connections (solid). Inter-slice connections (dashed) are encoded 
by Cj rs , specifying coupling of node j to itself between slices r and s. For clarity, inter- slice 
couplings are shown for only two nodes and depict two different types of couplings: (1) cou- 
pling between neighboring slices, appropriate for ordered slices; and (2) all-to-all inter-slice 
coupling, appropriate for categorical slices. 

Fig. 2: Multislice community detection of the Zachary Karate Club network (22) across mul- 
tiple resolutions. Colors depict community assignments of the 34 nodes (renumbered verti- 
cally to group similarly-assigned nodes) in each of the 16 slices (with resolution parameters 
7 S = {0.25, 0.5, . . . , 4}), for u = (top), u = 0.1 (middle), and u — 1 (bottom). Dashed lines 
bound the communities obtained using Newman-Girvan modularity (12). 

Fig. 3: Multislice community detection of U.S. Senate roll call vote similarities (23) with 
u = 0.5 coupling of 110 slices across time (110 two-year Congresses, covering 1789-2008). 
(A) Colors indicate assignments to 9 communities of the 1884 unique Senators (sorted verti- 
cally and connected across Congresses by dashed lines) in each Congress they appear. The dark 
blue and red communities correspond closely to the modern Democratic and Republican parties, 
respectively. Horizontal bars indicate the historical period of each community, with accompa- 
nying text enumerating nominal party affiliations of the single- slice nodes (each representing 
a Senator in a Congress): Pro-Administration (PA), Anti-Administration (AA), Federalist (F), 
Democratic-Republican (DR), Whig (W), Anti- Jackson (AJ), Adams (A), Jackson (J), Demo- 
cratic (D), and Republican (R). Vertical gray bars indicate Congresses in which three commu- 
nities appeared simultaneously. (B) The same assignments according to state affiliations. 
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We provide additional details here about the results and examples that we discussed in the 
main text. We begin by reviewing salient results from Ref. (13) concerning the connection 
between normalized Laplacian dynamics on networks and the modularity quality function for 
network community structure. We generalized this methodology to reproduce the null models 
for bipartite, directed, and signed networks, culminating in the specification of the correspond- 
ing quality function for multislice networks, which are combinations of individual networks 
coupled through links that connect each node in one network slice to itself in other slices. The 
stacking of multiple slices, linked together by identity arcs, provides a useful representation 
for visualization and extension of network measures to dynamic graphs (SI). We developed 
a methodology for community detection in such multislice networks, derived from stability of 
communities under normalized Laplacian dynamics. We additionally considered a similar anal- 
ysis following from the standard (unnormalized) Laplacian dynamics results from Ref. (13). 
We proved that the domains of optimization of each network partition are convex in the space 
of parameters for quality functions that are linear in those parameters and comment on possible 
consequences of this result. 

Laplacian Dynamics Formalism 

We review the Laplacian dynamics formalism recently developed by Lambiotte et al. (13). 

The crucial insight of Ref. (13) was to rederive network modularity from the continuous- 
time normalized Laplacian dynamics pi = J2j yM^j - ft on a unipartite, undirected network 
defined by the adjacency matrix components with node strengths fcj = J2j ^-ij- They also 
introduced a notion of stability of communities under such dynamics (13, S2) by directly com- 
paring the joint probability at stationarity of independent appearances at nodes % and j with 
the linear (in time) approximate map from node j to node i. In so doing, they derived a qual- 
ity function equivalent to Newman-Girvan (NG) modularity (12) at unit time and the standard 
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Potts generalization of NG modularity (5) that includes a resolution parameter (which is then 
interpreted as an inverse time). 

The normalized Laplacian dynamics, p\ = J2j ^^-ijPj ~ V ; n nave a steady state given by 
p* = kj/(2m), where 2m = J2ih = J2ij Mj describes the total strength (i.e., total edge weight) 
in the network. In Ref. (13), Lambiotte et al. quantified a measure of the stability R(t) of a 
specified partition of the network into communities using the probability that a random walker 
remains within the same community after time t, in statistically steady conditions, relative to 
that expected under independence. Using the operator Lij = Ay/kj — Sij of the dynamics, 
where 6% is the Kronecker delta, they specified this stability as 



R(t) = E 
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where the contribution from an independence assumption appears in the second term in brack- 
ets. Expanding the matrix exponential in equation ([!]) to first-order in t, so that (e* L )y 
5ij + tLij, Lambiotte et al. demonstrated that R(i) directly yields the quality function (13) 
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up to 5ij factors that always contribute to the sum and are thus immaterial in identifying parti- 
tions that optimize Q(t). The resulting quality function reduces to NG modularity for t = 1. 
Moreover, they showed that dividing by t (which has no effect on the optima for specified t) 
provides a direct interpretation of the resolution parameter 7 = 1/t when the quality is written 



in the usual form (5): Q = Y^ij 



A - — ryhihi 
^ ' 2m 



5(gi,gj). Hence, the stability of the commu- 



nity partition relative to that expected under independence provides a natural definition for the 
null model employed in the quality function. 
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Generalized Laplacian Dynamics 

We extended the formalism of Lambiotte et al. (13) to multislice networks by considering three 
crucial generalizations. 

First, we restricted the expected independent contribution given by the probability of a ran- 
dom walker remaining within the same community after time t in the statistically steady state 
to one that is conditional on the type of connection necessary to step between two nodes. That 
is, we replaced the p*p* independent contribution in equation (|TJ) with a conditional indepen- 
dent contribution Pi\jP*j, where pi\j is the conditional probability at stationarity of jumping to 
node i from node j along a specific edge type that is allowed by the specified category of net- 
works. This constraint on the independent contribution is consistent with the linear-in-time 
expansion of the exponential map employed in the calculation of the expected joint population, 
giving R(t) = Y,ij + tLij) p* — Pi\jP*j S{gi, 9j) after linearization of the exponential. For 
instance, we considered undirected bipartite networks. Such networks have two types of nodes 
(e.g., a person might belong to an organization), and every edge must connect a node of one 
type to a node of the other. The adjacency matrix A in the operator L^- = A+j/kj — Sij takes a 
bipartite form (i.e., it consists of off-diagonal blocks), and gives the same formula for the steady 
state, p* = kj/ (2m). However, = bijki/m, where is an indicator of bipartiteness that is 
equal to 1 if nodes i and j are of different types and otherwise. The denominator in is m 
(cf. 2m) because the probability of stepping to node i conditional on the additional information 
that the jump is along an edge going towards a node of z's type doubles the probability. Again 
neglecting 5^ contributions, dividing by t, and setting 7 = \/t, we thus obtained 



^bipartite 2 m X! 



m 



8(9*, 9 j), (3) 



which is the generalization of the (7 = 1) Barber bipartite null model (15) obtained by incor- 
porating the resolution parameter 7. 



Second, we generalized the Laplacian dynamics to include motion along multiple types of 
connections. For example, we considered a directed network (so that A is no longer symmetric) 
with kf' = J2j Aij and k° ut = J2i^-ij- We defined the normalized Laplacian dynamics to 
include motion equally along both incoming and outgoing edges, subject to the normalization 
kj = k™ 1 + kj Ut . That is, we studied the dynamics p\ = J2j LijPj = Y,j £r(Aj + ^ji)Pj ~ Pu 
which again has steady state p* = kj/(2m), with 2m = J2j kj = 2J2ij^-ij- The change 
induced by the consideration of the directed network occurs in the conditional probability p^, 
which must respect the type of edge (incoming versus outgoing) that is used to arrive at node i 
as well as the fraction of such edges available in the departure from node j. We thus obtained 



where each additive term combines the probability of picking a particular type of edge when 
departing node j with the probability of arriving at node % given that the motion is on that type 
of edge. Because of the symmetry in summing over {i, j} pairs, we have equivalently rewritten 
the resulting partition quality as 



which (as with bipartite networks) yielded the natural extension of the corresponding standard 
(7 = 1) null model for directed networks (16,17) by incorporation of the resolution parameter 7. 
This approach contrasts with that of Lambiotte et al. , which restricted consideration to motion 
following the link directions (13, S3). 

We also studied the Laplacian dynamics given by the operator = (Afj + A~j)/kj — 5ij 
(with kj — kj~ + kj), which similarly yielded a null model for (undirected) signed networks, 
in which link weights can be either positive or negative. Edges can be separated into ones 
that contribute positively (for which Afj > 0) and those that contribute negatively (for which 
Aj > 0). As before, the steady state is given by p* = kj/ (2m), where now m = m + + mr . 
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Because of the penalizing contribution desired from the AZ > links, we chose to weight 
the A~ and k~ contributions negatively when they appeared in the partition stability formula, 
which is given by equation ([TJ). Aside from this sign convention, we calculated the conditional 
probability at stationarity using the same procedure as in the directed case, keeping track of 
whether the movement from node j to node i is along a positive or negative edge. This gave 

* = ( kt kf _ K_kj_\ k^ = j L _ f ktkj _ k-k~ \ 

PlbPj * 2m+ kj 2m- kj I 2m 2m I 2m+ 2m- ' ' 



and yielded Q = ~ £ >3 
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S(gi,gj). This quality function re- 



duces at 7 = 1 to one proposed signed null model (18) and is a special case of a more general 
signed null model (19) that includes separate resolution parameters (j + and 7") for the positive 
and negative contributions. We reconstructed the latter null model by using our third general- 
ization, which we present next. 

Our third generalization was to flexibly interpret the stability under Laplacian dynamics in 
order to permit different spreading weights on the different types of links. This was not an issue 
in our consideration of directed networks unless one wants to weight incoming and outgoing 
edges differently. On the other hand, it might be desirable for signed networks to consider 
reweighted conditional probabilities at stationarity using some factor other than the relative 
strengths of the different edges at node j (even though we only considered a single specification 
of the underlying Laplacian dynamics). This generalization gave 
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S(gi,9j), (7) 



2m + 2m~ 

with two resolution parameters (7+ and 7""), which is the undirected version of the aforemen- 
tioned more general null model for signed networks (19), with the full directed version similarly 
obtained by combining the above generalizations. 

Having shown that our generalizations recovered the appropriate null models for other cat- 
egories of networks (bipartite, directed, and signed), we applied this methodology to the far 
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more general framework of multislice networks. We supposed that each slice s of a network is 
represented by adjacencies Aij S between nodes i and j and specified inter-slice couplings Cj rs 
that connect node j in slice r to itself in slice s (see Fig. 1). That is, notationally, we used two 
indices to specify each node-slice: a single node (e.g., i) in an indicated slice (e.g., s). For 
simplicity, we restricted our attention to undirected network slices (A^g = Aj is ) and undirected 
couplings (Cj rs = Cj sr ), but we can incorporate additional structure in the slices and couplings 
in the same manner as in the single-slice derivations above. For convenience, we notated the 
strengths of each node individually in each slice, so that kj S = Y,iAij S , Cj S = ^ r Cj sr , and 
we defined the multislice strength Kj S = k js + Cj S . We studied a continuous-time Laplacian 
process analogous to those above that respects the intra-slice nature of Aij S and the inter-slice 
couplings of Cj sr , specified by p is = J2j r (Aij S 5 S r + 5ijCj sr )pjr/ Kj r — p is , which has steady state 
probability distribution p* r = Kj r / (2/i), where 2/i = J2j r K jr- We then specified the associated 
multislice null model using the probability p is \j r of sampling node-slice (i, s) conditional on 
whether the multislice structure allows one to step from node-slice (j, r) to node-slice (i, s), 
considering intra- and inter-slice steps separately: 



The first term in brackets above describes the conditional probability appropriate for motion 
along intra-slice edges, analogous to those in the generalized derivation of other null models 
above, including the probability, kj r /Kj r , of using an intra-slice edge when leaving (j, r) and 
the resulting restriction to the given slice (5 sr ) made explicit. The second term in brackets, 
which similarly describes the conditional probability of motion between two slices, leverages 
the known definition of the Cj sr coupling. That is, the conditional probability of stepping 
from (j, r) to (i, s) along an inter-slice coupling is non-zero if and only if i = j, and it is 
proportional to the probability Cj sr /K,j r of selecting the precise inter-slice link that connects 



Pis\jrPj r 



kie kj 



® sr i 0. 
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to slice s from all edges connected to (j,r). The inter-slice strengths Cj r therefore canceled 
naturally as part of this calculation. Subtracting this conditional joint probability from the linear 
(in time) approximation of the exponential describing the Laplacian dynamics on the multislice 
networks, we then obtained a multislice generalization of modularity: 

Q multislice 7) ^ ' | I ^-ijs Ts - ^ J &sr &ijCj sr / S^Qisi 9jr) i (9) 

l ^ ijsr { \ lm s J ) 

where we have again utilized reweighting of the conditional probabilities, allowing for differ- 
ent resolutions 7 S in each slice. We absorbed the corresponding resolution parameter for the 
inter-slice couplings into the magnitude of the elements of Cj sr , which we then supposed for 
simplicity take binary values {0, cu} indicating absence/presence of inter-slice links. 

In the absence of such a reweighting in the interpretation of the stability of the partition, 
with 7s = 7 for all s, the corresponding prefactor on Cj sr absorbed above is (1 — 7). Impos- 
ing the choice 7 = 1 then recovered the usual interpretation of modularity as a count of the 
total weight of intra-slice edges minus the weight expected at random, and (as expected) the 
specified deterministic Cj sr contribution dropped out entirely, because such inter-slice links are 
definitional to the multislice network. In contrast, by leveraging the notion of stability under 
Laplacian dynamics, generalized appropriately, we have derived a principled generalization of 
modularity to multislice networks. 

Choosing binary-valued Cj sr = {0, cu} requires only a single coupling parameter u to con- 
trol the extent of inter-slice correspondence of communities. When to = 0, there is no benefit 
from extending communities across slices, so the optimal partition is obtained from independent 
optimization of the corresponding quality function in each slice. At the other extreme, when u 
becomes sufficiently large, the quality-optimizing partitions force the community assignment of 
a node to remain the same across all slices in which that node appears, and the multislice qual- 
ity reduces to a difference between the adjacency matrix summed over the contributions from 
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the individual slices and the sum over the separate single- slice null models (with selected 7 S ). 
That is, the null model obtained in the limit of large to is not the same as the standard NG null 
model on the adjacency matrix summed across slices, which only relies on the total summed 
degrees; rather, the required sum of the single-slice null models respects the degree sequences 
of these different contributions separately. The generality of this framework also allows one 
to consider different weights across the Cj sr couplings, if deemed appropriate for a particular 
application. Additionally, we note that the linearity of equation ([£]) with respect to the {7^, u} 
parameters necessitates that the modularity-optimizing domain of a single partition is convex in 
this parameter space (as derived below, with a brief discussion of consequences). 

Unnormalized Multislice Laplacian Dynamics 

As discussed by Lambiotte et al. (13), a similar analysis of the stability of communities under 
standard (i.e., unnormalized) Laplacian dynamics can be used to yield a quality function with 
a null model corresponding to a uniform random graph (5). We generalized this result to the 
multislice setting using a natural definition of the relevant independent probabilities subject to 
conditions imposed by the network structure specific to our multislice setting (similar to our 
derivation for normalized Laplacian dynamics). 

We specified the standard Laplacian dynamics on a multislice network defined by A^ s and 
C jsr by pis = T,jr( A ijsS sr + SijC jsr )p jr / (k) - Pi s K is /(K), where angled brackets denote an 
average over the entire multislice network and we recall that Kj S = kj s + Cj S is the multislice 
strength. The steady-state probability distribution under these dynamics is constant. Hence, 
p* r = l/N, where N is the total number of nodes summed across slices in the multislice 
network. We then scale the conditional probability p is \j r of stepping from node j at slice r to 
node i at slice s appropriate to the selected standard dynamics, where the rate of leaving node j 
at slice r is proportional to Kj r [cf . the constant rate of leaving (j, r) in the normalized Laplacian 



8 



dynamics in the rest of this paper]. Given that we repeated the procedure of allowing different 
resolution parameters (inverse times) both within and across slices, it was sufficient for us to 
consider the conditional independent probability in the form 

Pis\jrP* jr = l$sr + C jsr 5ij] — . (10) 

Ignoring 5ij8 sr contributions to quality, which have no effect on identifying the optimal parti- 
tion, we obtained 

Q = Y1 ~ 7s) $sr + SijCjsr} 5(g iS} g jr ) (1 1) 

ijsr 

as the multislice generalization of the uniform random null model. Note that we once again 
absorbed the inter-slice coupling strength directly into the binary values of Cj sr = {0,o;}. 
(Again, if desired, one can also consider different weights across the Cj sr couplings.) As with 
the multislice null model that we obtained from normalized Laplacian dynamics, the limiting 
behaviors of this quality function are towards independent partitioning of each slice as cu — > 
and towards averaging over slices for u 1, though the latter is greatly simplified here since 
it is merely a sum over constant contributions, in contrast with the more detailed null model in 
the large coupling limit corresponding to normalized Laplacian dynamics. 

Convex Domains of Optimization 

We proved that the linearity of equation ^ with respect to the {7,,, to} parameters necessitates 
that the quality-optimizing domain of a single partition be convex in this parameter space. This 
result holds more generally for any community-detection quality function that is linear in its 
parameters. That is, if an identified partition of the network is the highest-quality partition at 
two points in parameter space, then it necessarily gives the best partition along the entire line 
segment connecting those two points. 

The proof of this convexity result followed from the consideration of a line in parameter 
space that contains two distinct optima at different points. For the purposes of this proof, we 
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notated the parameters (e.g., resolution parameters and/or inter-slice coupling strengths) by the 
vector array A and the modularity-like quality function as 

Q = J2 B ij$(9i, 9j) = i A ij ~ \ p ijp\ S (9h 9j) , (12) 

ij ijp 

where the notation i and j for the node indices naturally generalizes over the complete multislice 
network. That is, Q = B : % = (A — A • P) : where x is me common-community indicator 
with elements 5(gi, gj) specific to the selected partition. The meaning of the double contractions 



(e.g., B : x) over indices and dot products over parameters (A • P) is clear from equation ( |T2[ ). 

We then assumed without loss of generality that a partition specified by Xi is me unique 
optimum for parameters Ai, with A\ — A : Xi an d Pi = P : Xi defined so that Q\ = 
Ai — Ai • Pi. If the distinct partition specified by X2 is strictly optimal to xi f° r parameters A 2 
(with analogous definitions for Q 2 = A 2 — X2 • P2), then we showed it must follow that 

Q 1 =A 1 -X 1 -P 1 >A 2 -\i- P 2 and Q 2 = A 2 - A 2 • P 2 > A 1 - A 2 • P x . (13) 

We combined these inequalities to yield (A 2 — Ai) • Pi > (A 2 — Ai) • P 2 . We then considered 
a vector array A 3 that is colinear with A 2 and Ai, so that A 3 = A 2 + /(A 2 — Ai) with / > 0, 
which yielded the result that the quality of the Xi and % 2 partitions at A 3 must satisfy 

A 2 -A 3 -P 2 = A 2 -A 2 -P 2 -/(A 2 -Ai)-P 2 > Ai-A 2 -Pi-/(A 2 -Ai)-Pi = Ai-A 3 -Pi . (14) 

That is, the partition X2 is necessarily of higher quality than Xi at A 3 (though neither of them 
needs to be the optimum there). Therefore, non-convex domains of optimization are forbidden 



in the parameter space of quality functions of the form in equation ( 12 ). 

This requirement of convex domains of quality optimization might be useful for comparing 
results across different resolution and coupling parameters, not only in the present multislice 
setting but for any network-partitioning quality function that is linear in resolution parameters. 
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Although other quality functions might of course be considered, we note that each quality func- 



tion discussed in the present manuscript is of the general form in equation ( 12 ). Computational 
results that do not conform to convex domains of optimization typically indicate regions in 
which further computation should uncover better optima. Indeed, for a particular application, 
it might be important to consider many different parameter choices in our generalized quality 
function. We do not worry about such details here, as our goal has been to present a framework 
that allows one to study the community structure of multislice networks, but it is neverthe- 
less important to mention it for further consideration. We additionally note that optimizing the 
standard modularity quality function is known to be an NP-complete problem (S4), and the 
cautionary observations regarding modularity optimization (20, 21) naturally also apply to our 
more general multislice framework. 

Examples 

We conclude by providing additional details for the three examples discussed in the main text. 

Community Detection Across Multiple Scales 

We performed simultaneous community detection across multiple resolutions (scales) in the 
well-known Zachary Karate Club benchmark network, which encodes the friendships between 
34 members of a karate club at a U.S. university in the 1970s (22). Keeping the same 34- 
node unweighted adjacency matrix across slices (so that Aij S = Ay for all s), the resolution 
associated with each slice is dictated by a value from a specified sequence of ^ s parameters, 
which we chose to be the 16 values 7 S = {0.25, 0.5, 0.75, . . . , 4}. In Fig. 2, we depict the 
community assignments that we obtained when the individual nodes are coupled with strengths 
to = {0, 0.1, 1} between each neighboring pair of the 16 ordered slices. For each u, we took the 
higher quality partition from that given by a spectral method plus Kernighan-Lin (KL) node- 
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swapping steps (4, 25) and a generalization of the Louvain algorithm (S5) plus KL steps. We 
note that, despite this approach, the depicted u — 1 partition can be clearly improved by lever- 
aging the definition of the inter-slice coupling; specifically, the communities of nodes 30-34 (in 
the renumbering in Fig. 2) at different resolutions can be merged to improve the total quality of 
the multislice partition. Future algorithmic improvements could explicitly identify similar situ- 
ations where merging or breaking communities across slices might improve the overall quality. 

When cj = 0, the optimal partition obtained corresponds to the union of the independent 
partitions of each separate resolution parameter. As u is increased, the coupling between neigh- 
boring slices encourages the partition to include communities that straddle multiple slices in 
the hierarchy of scales. The mathematical limit of arbitrarily large cu requires that, eventually, 
the communities span the full range of the considered resolutions. Because only the resolution 
parameters differed from one slice to the next in this multiple-resolution example, the limit of 
infinitely large inter-slice coupling here corresponded to single-resolution community detection 
at the average of the selected 7 S values, (7^) rs 2.125. Even at the smallest value of the resolu- 
tion parameter that we used (7 = 0.25), we already observed a split into two communities when 
00 > (recalling that the actual club fractured into two groups). We simultaneously obtained all 
of the other network scales, such as the partitioning of the Karate Club into four communities at 
the default resolution of NG modularity (3, 25). We also identified nodes that have an especially 
strong tendency to break off from larger communities (e.g., nodes 24-29 in Fig. 2). 

This example illustrated that multislice community detection makes it possible to systemat- 
ically track the development of multiple network scales simultaneously. 

Community Detection in Time-Dependent Networks 

We considered roll call voting in the United States Senate across time. The Senate is one of 
the two chambers of the legislative branch (collectively called the Congress) of the U.S. federal 
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government. It currently consists of 100 Senators (two from each state) who serve staggered 
six-year terms such that approximately one-third of the Senate is elected every two years. The 
data we studied is from the lst-HOth Congresses, covering the years 1789-2008 and includ- 
ing 1884 individual Senators^ With each slice (i.e., within each two-year Congress), we de- 
fined a weighted connection between each pair of Senators in terms of a similarity between the 
votes they cast during that Congress (23). We then demonstrated that one can gain additional 
understanding of this network, and the underlying political processes, by applying multislice 
community detection to the collection of these 110 network slices taken as a whole. In this 
multislice network, we coupled each individual Senator to him/herself when appearing in con- 
secutive Congresses. If a Senator from Congress s did not serve in Congress s + 1, then we 
did not introduce inter-slice coupling between slices s and s + 1 for this individual. With this 
formulation, link strengths and nodes (Senators) both changed from one slice to another. 

Multislice community detection uncovered details about the individual and group voting 
dynamics over time that are simply not captured by the union of the 110 independent partitions 
of the individual Congresses. Again using a generalization of the Louvain algorithm plus KL 
steps, and using inter-slice coupling tu = 0.5, we obtained the partition depicted in Fig. 3 of 
the 1884 unique U.S. Senators in each Congress in which they voted into 9 communities. This 
community structure highlights several historical turning points in U.S. politics. For instance, 
the Congresses in which three communities appeared simultaneously are each historically sig- 
nificant: The 4th and 5th Congresses were the first with political parties; the 10th and 11th 
Congresses occurred during the political drama of former Vice President Aaron Burr's indict- 
ment for treason; the 14th and 15th Congresses witnessed the beginning of changing group 
structures in the Democratic-Republican party (23) amidst the dying Federalist party; the 31st 

'At least five Senators in the data [available at voteview . com| (26)] are each assigned two different identi- 
fication numbers, corresponding to different periods of their careers. We take the data as provided, counting such 
Senators twice, and merely remark that politically-minded studies should include such considerations. 



13 



Congress included the Compromise of 1850; the 37th Congress occurred during the beginning 
of the American Civil War; the 73rd and 74th Congresses followed the landslide 1932 elec- 
tion amidst the Great Depression; and the 85th-88th Congresses brought the major American 
civil rights acts, including the Congressional fights over the Civil Rights Acts of 1957, 1960, 
and 1964 (observe that all 44 Democratic Senators in the community colored green during 
this time period came from Southern states). A more complete political study using multislice 
community detection, which includes systematically examining the community structure as the 
inter-slice coupling strength lu is varied, would enable one to investigate such observations in 
extensive detail. 

Community Detection in Multiplex Networks 

We applied multislice community detection to a multiplex network of 1640 college students at 
an anonymous, northeastern American university (24). We included the symmetrized connec- 
tions from the first wave of this data (covering the first year of university attendance) repre- 
senting (1) Facebook friendships; (2) picture friendships, in which a student posted and tagged 
a photograph of another online; (3) roommates, in which two students shared a first-year dor- 
mitory room, creating clusters of 1-6 students; and (4) "housing group" preferences identified 
by the students. Because the different tie types are categorical, the natural inter-slice couplings 
connect an individual corresponding to one type of connection to him/herself in each of the 
other 3 types of networks. This type of inter-slice coupling thus has a different nature from the 
inter-slice couplings above that connected only neighboring (ordered) network slices. 

In Table [TJ we provide a summary of the basic results that we obtained by varying the inter- 
slice coupling strength u>. We tabulated the total number of communities and the percentages 
of individuals assigned to 1,2, 3, or 4 communities in the multislice network across the four 
types of connections. Again, to = yielded separate communities for each slice, as expected, 
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with each individual placed into four separate communities. As to was increased, communities 
merged across slices — most predominantly where the patterns of connection were relatively 
similar between two slices. This reduced the total number of communities and resulted in indi- 
viduals with fewer distinct community assignments across their 4 appearances in the different 
network slices. For u G [0.2,0.5], a significant majority of the individuals were assigned to 
only 1 or 2 communities, indicating that their social networks maintain group-level similarities 
across the four types of connections. Another significant set of students were grouped into 3 dif- 
ferent communities, and a small minority maintained 4 separate assignments, suggesting stark 
differences in their positions in the 4 single-category network slices. Finally, for u — 1, the 
inter-slice coupling was sufficiently strong that it forced all 4 multislice nodes corresponding to 
an individual student to be assigned to the same community. Further investigation of such dif- 
ferent community assignments across slices could be used to more clearly compare and contrast 
the roles of individuals in each network slice and in the complete multislice network. Addition- 
ally, a multislice approach might provide a novel mechanism for dealing with the problem of 
overlapping community assignments, as the hard partitioning of each node (located in a single 
slice) in the multislice network allows an individual to be placed into different communities in 
their appearances in the different slices. Indeed, multiplexity is itself a strong motivation for 
developing methods that allow communities to overlap (2,3,S6). 
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